	program xelf
c     test driver for "elf"
	real*8 PI,PIHALF
	parameter (PI=3.1415926535897932384626433d0)
	parameter (PIHALF=1.5707963267948966192313216916398d0)
	real*8 dmc,dphi,mc,mm,phi,phic,xf,elf
	real*4 rmc,rphi,rphic,rf,relf
	integer jend,iend,j,i
c
	jend=6
	iend=5
	dmc=1.d0/dble(jend-1)
	dphi=PIHALF/dble(iend)
 	write(*,'(1x,2a10,a25,a15)')	'm','phi/PI','elf','relf'
	do j=1,jend
		write(*,'(1x)')
		mc=dble(j-1)*dmc
		if(mc.le.0.d0) mc=1.21d-32
		rmc=mc
		mm=1.d0-mc
		do i=0,iend
		  phi=dphi*dble(i)
		  phic=dphi*dble(iend-i)
		  rphi=phi
	    rphic=phic
		  xf=elf(phi,phic,mc)
		  rf=relf(rphi,rphic,rmc)
          write(*,'(1x,0p2f10.5,0p1f25.15,0p1f15.7)')
     &      mm,phi/PI,xf,rf
		enddo
	enddo
	END
c---------------------------------------------------------------------------
      real*8 function elf(phi,phic,mc)
c
c	Double precision incomplete elliptic integral of the first kind
c
c     Reference: T. Fukushima, (2010) Numer. Math. 116, 687-719
c        "Fast Computation of Incomplete Elliptic Integral of First Kind
c         by Half Argument Transformation"		
c
c     Author: T. Fukushima Toshio.Fukushima@nao.ac.jp
c
c     Used subprograms: asn, acn, elk, serf (called from asn)
c
c     Inputs: phi  = argument                0 <= phi  <= PI/2
c             phic = complementar argument   0 <= phic <= PI/2
c             mc   = complementary parameter 0 <= mc   <= 1
c
c     Output: elf
c
c     CAUTION: phi and phic must satisfy condition, phi + phic = PI/2
c
      real*8 phi,phic,mc
	real*8 asn,elk,acn
      real*8 m,c,yc,d2,v
c
	m=1.d0-mc
	if(phi.lt.1.25d0) then
		elf=asn(sin(phi),m)
	else
		c=sin(phic)
		yc=c*c
		d2=mc+m*yc
		if(yc.lt.0.9d0*d2) then
			elf=elk(mc)-asn(c/sqrt(d2),m)
		else
			v=mc*(1.d0-yc)
			if(v.lt.yc*d2) then
				elf=acn(c,mc)
			else
				elf=elk(mc)-asn(sqrt(v/d2),m)
			endif
		endif
	endif
	return
      end
c---------------------------------------------------------------------------
      real*8 function elk(mc)
c
c	Double precision complete elliptic integral of the first kind
c
c     Reference: T. Fukushima, (2009) Celest. Mech. Dyn. Astron. 105, 305-328
c        "Fast Computation of Complete Elliptic Integrlals and Jacobian
c         Elliptic Functions"
c
c     Author: T. Fukushima Toshio.Fukushima@nao.ac.jp
c
c     Inputs: mc   = complementary parameter 0 <= mc   <= 1
c
c     Output: elk
c
	real*8 mc
	real*8 mcold,PIHALF,PIINV,elkold,TINY,m,mx,P,Q
	real*8 kkc,nome
c
	real*8 D1,D2,D3,D4,D5,D6,D7,D8,D9,D10,D11,D12,D13,D14
	parameter (D1=1.d0/16.d0,D2=1.d0/32.d0,D3=21.d0/1024.d0)
	parameter (D4=31.d0/2048.d0,D5=6257.d0/524288.d0)
	parameter (D6=10293.d0/1048576.d0,D7=279025.d0/33554432.d0)
	parameter (D8=483127.d0/67108864.d0)
	parameter (D9=435506703.d0/68719476736.d0)
	parameter (D10=776957575.d0/137438953472.d0)
	parameter (D11=22417045555.d0/4398046511104.d0)
	parameter (D12=40784671953.d0/8796093022208.d0)
	parameter (D13=9569130097211.d0/2251799813685248.d0)
	parameter (D14=17652604545791.d0/4503599627370496.d0)
c
	logical first/.TRUE./
	save first,mcold,PIHALF,PIINV,elkold,TINY
c
	if(first) then
		first=.FALSE.
		mcold=1.d0
		PIHALF=atan(1.d0)*2.d0
		PIINV=0.5d0/PIHALF
		elkold=PIHALF
		TINY=1.d-99
	endif
	m=1.d0-mc
	if(abs(m).lt.1.d-16) then
		elk=PIHALF
	elseif(abs(mc-mcold).lt.1.11d-16*mc) then
		elk=elkold
	elseif(mc.lt.TINY) then
	  elk=1.3862943611198906d0-0.5d0*log(TINY)
	elseif(mc.lt.1.11d-16) then
	  elk=1.3862943611198906d0-0.5d0*log(mc)
	elseif(mc.lt.0.1d0) then
		nome=mc*(D1+mc*(D2+mc*(D3+mc*(D4+mc*(D5+mc*(D6
     &		+mc*(D7+mc*(D8+mc*(D9+mc*(D10+mc*(D11+mc*(D12
     &		+mc*(D13+mc*D14)))))))))))))
		mx=mc-0.05d0
c
c	K'
c
		kkc=1.591003453790792180d0+mx*(
     &		0.416000743991786912d0+mx*(
     &		0.245791514264103415d0+mx*(
     &		0.179481482914906162d0+mx*(
     &		0.144556057087555150d0+mx*(
     &		0.123200993312427711d0+mx*(
     &		0.108938811574293531d0+mx*(
     &		0.098853409871592910d0+mx*(
     &		0.091439629201749751d0+mx*(
     &		0.085842591595413900d0+mx*(
     &		0.081541118718303215d0))))))))))
c
		elk=-kkc*PIINV*log(nome)
	elseif(m.le.0.1d0) then
		mx=m-0.05d0
		elk=1.591003453790792180d0+mx*(
     &		0.416000743991786912d0+mx*(
     &		0.245791514264103415d0+mx*(
     &		0.179481482914906162d0+mx*(
     &		0.144556057087555150d0+mx*(
     &		0.123200993312427711d0+mx*(
     &		0.108938811574293531d0+mx*(
     &		0.098853409871592910d0+mx*(
     &		0.091439629201749751d0+mx*(
     &		0.085842591595413900d0+mx*(
     &		0.081541118718303215d0))))))))))
	elseif(m.le.0.2d0) then
		mx=m-0.15d0
		elk=1.635256732264579992d0+mx*(
     &		0.471190626148732291d0+mx*(
     &		0.309728410831499587d0+mx*(
     &		0.252208311773135699d0+mx*(
     &		0.226725623219684650d0+mx*(
     &		0.215774446729585976d0+mx*(
     &		0.213108771877348910d0+mx*(
     &		0.216029124605188282d0+mx*(
     &		0.223255831633057896d0+mx*(
     &		0.234180501294209925d0+mx*(
     &		0.248557682972264071d0+mx*(
     &		0.266363809892617521d0+mx*(
     &		0.287728452156114668d0))))))))))))
	elseif(m.le.0.3d0) then
		mx=m-0.25d0
		elk=1.685750354812596043d0+mx*(
     &		0.541731848613280329d0+mx*(
     &		0.401524438390690257d0+mx*(
     &		0.369642473420889090d0+mx*(
     &		0.376060715354583645d0+mx*(
     &		0.405235887085125919d0+mx*(
     &		0.453294381753999079d0+mx*(
     &		0.520518947651184205d0+mx*(
     &		0.609426039204995055d0+mx*(
     &		0.724263522282908870d0+mx*(
     &		0.871013847709812357d0+mx*(
     &		1.057652872753547036d0)))))))))))
	elseif(m.le.0.4d0) then
		mx=m-0.35d0
		elk=1.744350597225613243d0+mx*(
     &		0.634864275371935304d0+mx*(
     &		0.539842564164445538d0+mx*(
     &		0.571892705193787391d0+mx*(
     &		0.670295136265406100d0+mx*(
     &		0.832586590010977199d0+mx*(
     &		1.073857448247933265d0+mx*(
     &		1.422091460675497751d0+mx*(
     &		1.920387183402304829d0+mx*(
     &		2.632552548331654201d0+mx*(
     &		3.652109747319039160d0+mx*(
     &		5.115867135558865806d0+mx*(
     &		7.224080007363877411d0))))))))))))
	elseif(m.le.0.5d0) then
		mx=m-0.45d0
		elk=1.813883936816982644d0+mx*(
     &		0.763163245700557246d0+mx*(
     &		0.761928605321595831d0+mx*(
     &		0.951074653668427927d0+mx*(
     &		1.315180671703161215d0+mx*(
     &		1.928560693477410941d0+mx*(
     &		2.937509342531378755d0+mx*(
     &		4.594894405442878062d0+mx*(
     &		7.330071221881720772d0+mx*(
     &		11.87151259742530180d0+mx*(
     &		19.45851374822937738d0+mx*(
     &		32.20638657246426863d0+mx*(
     &		53.73749198700554656d0+mx*(
     &		90.27388602940998849d0)))))))))))))
	elseif(m.le.0.6d0) then
		mx=m-0.55d0
		elk=1.898924910271553526d0+mx*(
     &		0.950521794618244435d0+mx*(
     &		1.151077589959015808d0+mx*(
     &		1.750239106986300540d0+mx*(
     &		2.952676812636875180d0+mx*(
     &		5.285800396121450889d0+mx*(
     &		9.832485716659979747d0+mx*(
     &		18.78714868327559562d0+mx*(
     &		36.61468615273698145d0+mx*(
     &		72.45292395127771801d0+mx*(
     &		145.1079577347069102d0+mx*(
     &		293.4786396308497026d0+mx*(
     &		598.3851815055010179d0+mx*(
     &		1228.420013075863451d0+mx*(
     &		2536.529755382764488d0))))))))))))))
	elseif(m.le.0.7d0) then
		mx=m-0.65d0
		elk=2.007598398424376302d0+mx*(
     &		1.248457231212347337d0+mx*(
     &		1.926234657076479729d0+mx*(
     &		3.751289640087587680d0+mx*(
     &		8.119944554932045802d0+mx*(
     &		18.66572130873555361d0+mx*(
     &		44.60392484291437063d0+mx*(
     &		109.5092054309498377d0+mx*(
     &		274.2779548232413480d0+mx*(
     &		697.5598008606326163d0+mx*(
     &		1795.716014500247129d0+mx*(
     &		4668.381716790389910d0+mx*(
     &		12235.76246813664335d0+mx*(
     &		32290.17809718320818d0+mx*(
     &		85713.07608195964685d0+mx*(
     &		228672.1890493117096d0+mx*(
     &		612757.2711915852774d0))))))))))))))))
	elseif(m.le.0.8d0) then
		mx=m-0.75d0
		elk=2.156515647499643235d0+mx*(
     &		1.791805641849463243d0+mx*(
     &		3.826751287465713147d0+mx*(
     &		10.38672468363797208d0+mx*(
     &		31.40331405468070290d0+mx*(
     &		100.9237039498695416d0+mx*(
     &		337.3268282632272897d0+mx*(
     &		1158.707930567827917d0+mx*(
     &		4060.990742193632092d0+mx*(
     &		14454.00184034344795d0+mx*(
     &		52076.66107599404803d0+mx*(
     &		189493.6591462156887d0+mx*(
     &		695184.5762413896145d0+mx*(
     &		2.567994048255284686d6+mx*(
     &		9.541921966748386322d6+mx*(
     &		3.563492744218076174d7+mx*(
     &		1.336692984612040871d8+mx*(
     &		5.033521866866284541d8+mx*(
     &		1.901975729538660119d9+mx*(
     &		7.208915015330103756d9)))))))))))))))))))
	elseif(m.le.0.85d0) then
		mx=m-0.825d0
		elk=2.318122621712510589d0+mx*(
     &		2.616920150291232841d0+mx*(
     &		7.897935075731355823d0+mx*(
     &		30.50239715446672327d0+mx*(
     &		131.4869365523528456d0+mx*(
     &		602.9847637356491617d0+mx*(
     &		2877.024617809972641d0+mx*(
     &		14110.51991915180325d0+mx*(
     &		70621.44088156540229d0+mx*(
     &		358977.2665825309926d0+mx*(
     &		1.847238263723971684d6+mx*(
     &		9.600515416049214109d6+mx*(
     &		5.030767708502366879d7+mx*(
     &		2.654441886527127967d8+mx*(
     &		1.408862325028702687d9+mx*(
     &		7.515687935373774627d9)))))))))))))))
      else
		mx=m-0.875d0
		elk=2.473596173751343912d0+mx*(
     &		3.727624244118099310d0+mx*(
     &		15.60739303554930496d0+mx*(
     &		84.12850842805887747d0+mx*(
     &		506.9818197040613935d0+mx*(
     &		3252.277058145123644d0+mx*(
     &		21713.24241957434256d0+mx*(
     &		149037.0451890932766d0+mx*(
     &		1.043999331089990839d6+mx*(
     &		7.427974817042038995d6+mx*(
     &		5.350383967558661151d7+mx*(
     &		3.892498869948708474d8+mx*(
     &		2.855288351100810619d9+mx*(
     &		2.109007703876684053d10+mx*(
     &		1.566998339477902014d11+mx*(
     &		1.170222242422439893d12+mx*(
     &		8.777948323668937971d12+mx*(
     &		6.610124275248495041d13+mx*(
     &		4.994880537133887989d14+mx*(
     &		3.785974339724029920d15)))))))))))))))))))
	endif
c
	mcold=mc
	elkold=elk
      return
      end
c---------------------------------------------------------------------------
      real*8 function acn(c0,mc)
c
      real*8 c0,mc
	real*8 m,c,yc,y,s,f,d
	real*8 asn
      integer i
	integer jsn,jcn
	common /count/ jsn,jcn
c
	m=1.d0-mc
	c=c0
	yc=c*c
	if(yc.gt.0.5d0) then
		y=1.d0-yc
		s=sqrt(y)
		acn=asn(s,m)
	jcn=0
		return
	endif
	f=1.d0
      do i=1,20
    		d=sqrt(mc+m*yc)
		yc=(c+d)/(1.d0+d)
		f=f*2.d0
		if(yc.gt.0.5d0) then
			y=1.d0-yc
			s=sqrt(y)
	jcn=i
			goto 1
		endif
		c=sqrt(yc)
	enddo
	write(*,*) "(acn) too many iterations: c0,mc=",c0,mc
    1 continue
	acn=f*asn(s,m)
	return
      end
c---------------------------------------------------------------------------
      real*8 function asn(s0,m)
c
      real*8 s0,m
      real*8 del,s,f,y
	real*8 serf
	integer j
	integer jsn,jcn
	common /count/ jsn,jcn
c
	del=0.04094d0-0.00652d0*m
	s=s0
	y=s*s
	if(y.lt.del) then
		asn=s*serf(y,m)
	jsn=0
		return
	endif
	f=1.d0
      do j=1,20
    		y=y/((1.d0+sqrt(1.d0-y))*(1.d0+sqrt(1.d0-m*y)))
		f=f*2.d0
		if(y.lt.del) then
			s=sqrt(y)
	jsn=j
			goto 1
		endif
	enddo
	write(*,*) "(asn) too many iterations: s0,m=",s0,m
    1 continue
	asn=f*s*serf(y,m)
	return
      end
c---------------------------------------------------------------------------
      real*8 function serf(y,m)
c
      real*8 y,m
c
	real*8 F1,F2,F3,F4
	real*8 F10,F20,F21,F30,F31,F40,F41,F42
	real*8 F5,F50,F51,F52,F6,F60,F61,F62,F63
	real*8 F7,F70,F71,F72,F73,F8,F80,F81,F82,F83,F84
	real*8 F9,F90,F91,F92,F93,F94
	parameter (F10=1.d0/6.d0)
	parameter (F20=3.d0/40.d0)
	parameter (F21=2.d0/40.d0)
	parameter (F30=5.d0/112.d0)
	parameter (F31=3.d0/112.d0)
	parameter (F40=35.d0/1152.d0)
	parameter (F41=20.d0/1152.d0)
	parameter (F42=18.d0/1152.d0)
	parameter (F50=63.d0/2816.d0)
	parameter (F51=35.d0/2816.d0)
	parameter (F52=30.d0/2816.d0)
	parameter (F60=231.d0/13312.d0)
	parameter (F61=126.d0/13312.d0)
	parameter (F62=105.d0/13312.d0)
	parameter (F63=100.d0/13312.d0)
	parameter (F70=429.d0/30720.d0)
	parameter (F71=231.d0/30720.d0)
	parameter (F72=189.d0/30720.d0)
	parameter (F73=175.d0/30720.d0)
	parameter (F80=6435.d0/557056.d0)
	parameter (F81=3432.d0/557056.d0)
	parameter (F82=2722.d0/557056.d0)
	parameter (F83=2520.d0/557056.d0)
	parameter (F84=2450.d0/557056.d0)
	parameter (F90=12155.d0/1245184.d0)
	parameter (F91=6435.d0/1245184.d0)
	parameter (F92=5148.d0/1245184.d0)
	parameter (F93=4620.d0/1245184.d0)
	parameter (F94=4410.d0/1245184.d0)
c
	F1=F10+m*F10
	F2=F20+m*(F21+m*F20)
	F3=F30+m*(F31+m*(F31+m*F30))
	F4=F40+m*(F41+m*(F42+m*(F41+m*F40)))
	F5=F50+m*(F51+m*(F52+m*(F52+m*(F51+m*F50))))
	F6=F60+m*(F61+m*(F62+m*(F63+m*(F62+m*(F61+m*F60)))))
	F7=F70+m*(F71+m*(F72+m*(F73+m*(F73+m*(F72+m*(F71+m*F70))))))
	F8=F80+m*(F81+m*(F82+m*(F83+m*(F84
     &	+m*(F83+m*(F82+m*(F81+m*F80)))))))
	F9=F90+m*(F91+m*(F92+m*(F93+m*(F94+m*(F94
     &	+m*(F93+m*(F92+m*(F91+m*F90))))))))
	serf=1.d0+y*(F1+y*(F2+y*(F3+y*(F4
     &	+y*(F5+y*(F6+y*(F7+y*(F8+y*F9))))))))
	return
      end
c---------------------------------------------------------------------------
      real function relf(phi,phic,mc)
c
c	Single precision incomplete elliptic integral of the fist kind
c
c     Reference: T. Fukushima, (2010) Numer. Math. 116, 687-719
c        "Fast Computation of Incomplete Elliptic Integral of First Kind
c         by Half Argument Transformation"		
c
c     Author: T. Fukushima Toshio.Fukushima@nao.ac.jp
c
c     Used subprograms: rasn, racn, relk, rserf (called from rasn)
c
c     Inputs: phi  = argument                0 <= phi  <= PI/2
c             phic = complementar argument   0 <= phic <= PI/2
c             mc   = complementary parameter 0 <= mc   <= 1
c
c     Output: relf
c
c     CAUTION: phi and phic must satisfy condition, phi + phic = PI/2
c
      real phi,phic,mc
	real rasn,relk,racn
      real m,c,yc,d2,v
c
	m=1.0-mc
	if(phi.lt.1.25) then
		relf=rasn(sin(phi),m)
	else
		c=sin(phic)
		yc=c*c
		d2=mc+m*yc
		if(yc.lt.0.9*d2) then
			relf=relk(mc)-rasn(c/sqrt(d2),m)
		else
			v=mc*(1.0-yc)
			if(v.lt.yc*d2) then
				relf=racn(c,mc)
			else
				relf=relk(mc)-rasn(sqrt(v/d2),m)
			endif
		endif
	endif
	return
      end
c---------------------------------------------------------------------------
      real function relk(mc)
c
c	Single precision complete elliptic integral of the first kind
c
c     Reference: T. Fukushima, (2009) Celest. Mech. Dyn. Astron. 105, 305-328
c        "Fast Computation of Complete Elliptic Integrlals and Jacobian
c         Elliptic Functions"
c
c     Author: T. Fukushima Toshio.Fukushima@nao.ac.jp
c
c     Inputs: mc   = complementary parameter 0 <= mc   <= 1
c
c     Output: elk
c
	real mc
	real m,nome,mx,kkc,P,Q
c
	real D1,D2,D3,D4,D5,D6
	parameter (D1=1.0/16.0,D2=1.0/32.0,D3=21.0/1024.0)
	parameter (D4=31.0/2048.0,D5=6257.0/524288.0)
	parameter (D6=10293.0/1048576.0)
	real PIHALF,PIINV
	parameter (PIHALF=1.57079633,PIINV=0.318309886)
c
	m=1.0-mc
	if(mc.lt.1.05e-8) then
	    relk=1.38629436-0.5*log(mc)
	elseif(mc.lt.0.1) then
		nome=mc*(D1+mc*(D2+mc*(D3+mc*(D4+mc*(D5+mc*D6)))))
		mx=mc-0.05
c
c	K'
c
		kkc=1.59100345+mx*(
     &		0.41600074+mx*(
     &		0.24579151+mx*(
     &		0.17948148+mx*(
     &		0.14455606))))
c
		relk=-kkc*PIINV*log(nome)
	elseif(m.le.0.1) then
		mx=m-0.05
		relk=1.59100345+mx*(
     &		0.41600074+mx*(
     &		0.24579151+mx*(
     &		0.17948148+mx*(
     &		0.14455606))))
	elseif(m.le.0.2) then
		mx=m-0.15
		relk=1.63525673+mx*(
     &		0.47119063+mx*(
     &		0.30972841+mx*(
     &		0.25220831+mx*(
     &		0.22672562))))
	elseif(m.le.0.3) then
		mx=m-0.25
		relk=1.68575035+mx*(
     &		0.54173185+mx*(
     &		0.40152444+mx*(
     &		0.36964247+mx*(
     &		0.37606072))))
	elseif(m.le.0.4) then
		mx=m-0.35
		relk=1.74435060+mx*(
     &		0.63486428+mx*(
     &		0.53984256+mx*(
     &		0.57189271+mx*(
     &		0.67029514+mx*(
     &		0.83258659)))))
	elseif(m.le.0.5) then
		mx=m-0.45
		relk=1.81388394+mx*(
     &		0.76316325+mx*(
     &		0.76192861+mx*(
     &		0.95107465+mx*(
     &		1.31518067+mx*(
     &		1.92856069)))))
	elseif(m.le.0.6) then
		mx=m-0.55
		relk=1.89892491+mx*(
     &		0.95052179+mx*(
     &		1.15107759+mx*(
     &		1.75023911+mx*(
     &		2.95267681+mx*(
     &		5.28580040)))))
	elseif(m.le.0.7) then
		mx=m-0.65
		relk=2.00759840+mx*(
     &		1.24845723+mx*(
     &		1.92623466+mx*(
     &		3.75128964+mx*(
     &		8.11994455+mx*(
     &		18.6657213+mx*(
     &		44.6039248))))))
	elseif(m.le.0.8) then
		mx=m-0.75
		relk=2.15651565+mx*(
     &		1.79180564+mx*(
     &		3.82675129+mx*(
     &		10.3867247+mx*(
     &		31.4033141+mx*(
     &		100.923704+mx*(
     &		337.326828+mx*(
     &		1158.70793)))))))
	elseif(m.le.0.85) then
		mx=m-0.825
		relk=2.31812262+mx*(
     &		2.61692015+mx*(
     &		7.89793508+mx*(
     &		30.5023972+mx*(
     &		131.486937+mx*(
     &		602.984764+mx*(
     &		2877.02462))))))
      else
		mx=m-0.875
		relk=2.47359617+mx*(
     &		3.72762424+mx*(
     &		15.6073930+mx*(
     &		84.1285084+mx*(
     &		506.981820+mx*(
     &		3252.27706+mx*(
     &		21713.2424+mx*(
     &		149037.045)))))))
	endif
c
      return
      end
c---------------------------------------------------------------------------
      real function racn(c0,mc)
c
      real c0,mc
	real m,c,yc,y,s,f,d
	real rasn
      integer i
	integer jsn,jcn
	common /count/ jsn,jcn
c
	m=1.0-mc
	c=c0
	yc=c*c
	if(yc.gt.0.5) then
		y=1.0-yc
		s=sqrt(y)
		racn=rasn(s,m)
		jcn=0
		return
	endif
	f=1.0
      do i=1,10
    		d=sqrt(mc+m*yc)
		yc=(c+d)/(1.0+d)
		f=f*2.0
		if(yc.gt.0.5) then
			y=1.0-yc
			s=sqrt(y)
			jcn=i
			goto 1
		endif
		c=sqrt(yc)
	enddo
	write(*,*) "(racn) too many iterations: c0,mc=",c0,mc
    1 continue
	racn=f*rasn(s,m)
	return
      end
c---------------------------------------------------------------------------
      real function rasn(s0,m)
c
      real s0,m
      real del,s,f,y
	real rserf
	integer j
	integer jsn,jcn
	common /count/ jsn,jcn
c
	del=0.1888-0.0378*m
	s=s0
	y=s*s
	if(y.lt.del) then
		rasn=s*rserf(y,m)
		jsn=0
		return
	endif
	f=1.0
      do j=1,10
    		y=y/((1.0+sqrt(1.0-y))*(1.0+sqrt(1.0-m*y)))
		f=f*2.0
		if(y.lt.del) then
			s=sqrt(y)
			jsn=j
			goto 1
		endif
	enddo
	write(*,*) "(rasn) too many iterations: s0,m=",s0,m
    1 continue
	rasn=f*s*rserf(y,m)
	return
      end
c---------------------------------------------------------------------------
      real function rserf(y,m)
c
      real y,m
c
	real F1,F2,F3,F4
	real F10,F20,F21,F30,F31,F40,F41,F42
	real F5,F50,F51,F52,F6,F60,F61,F62,F63
	parameter (F10=1.0/6.0)
	parameter (F20=3.0/40.0)
	parameter (F21=2.0/40.0)
	parameter (F30=5.0/112.0)
	parameter (F31=3.0/112.0)
	parameter (F40=35.0/1152.0)
	parameter (F41=20.0/1152.0)
	parameter (F42=18.0/1152.0)
	parameter (F50=63.0/2816.0)
	parameter (F51=35.0/2816.0)
	parameter (F52=30.0/2816.0)
	parameter (F60=231.0/13312.0)
	parameter (F61=126.0/13312.0)
	parameter (F62=105.0/13312.0)
	parameter (F63=100.0/13312.0)
c
	F1=F10+m*F10
	F2=F20+m*(F21+m*F20)
	F3=F30+m*(F31+m*(F31+m*F30))
	F4=F40+m*(F41+m*(F42+m*(F41+m*F40)))
	F5=F50+m*(F51+m*(F52+m*(F52+m*(F51+m*F50))))
	F6=F60+m*(F61+m*(F62+m*(F63+m*(F62+m*(F61+m*F60)))))
	rserf=1.0+y*(F1+y*(F2+y*(F3+y*(F4
     &	+y*(F5+y*F6)))))
	return
      end
      